Abnormality finding in projection images

ABSTRACT

The present invention includes a utility for determining the severity of a stenosis in a blood vessel. In one aspect, a method for improving DSA image quality includes: (1) registration of the mask and bolus images prior to subtracting procedure to reduce the artifacts from misalignment; (2) enhancement of the registered DSA image by an anisotropic diffusion technique and a nonlinear normalization technique; and (3) detecting the boundary of blood vessel and quantitatively measuring percentage stenosis, which may be done automatically.

CROSS-REFERENCE TO RELATED APPLICATION

This application claims priority under 35 U.S.C. §119 to U.S.Provisional Application No. 61/057,725 having a filing date of May 30,2008, the entire contents of which are incorporated by reference herein.

FIELD

The present disclosure is directed to medical imaging systems. Morespecifically, the present disclosure is directed to systems and methodsthat alone or collectively facilitate real-time imaging.

BACKGROUND

Coronary artery disease causes in excess of 1.5 million cases ofmyocardial infarction annually, and is the leading cause of death in theUnited States, resulting in more than 500,000 deaths per year. Theaccurate diagnosis and quantification of coronary artery disease iscritical to subsequent treatment decisions. Despite the emergence of newtechniques for the visualization and analysis of vascular structures,digital subtraction angiography (DSA) remains the preferred procedure tohelp clinicians to make clinical-decisions in all kinds of vasculardiseases. Digital subtraction angiography (DSA) is a well-establishedmodality for the visualization of blood vessels in the human body.

During a DSA procedure a first series of frames are taken whichrepresents the body anatomy in static form. They are called mask frames.Then dye is injected in the body through a catheter to visualize theblood vessel. The dye-injected frames are called the bolus frames. Asubtraction procedure is preformed between the averaged bolus frame andthe averaged mask frame to get DSA image. The patient is positioned onthe X-ray imaging system while an X-ray movie is acquired and a DSAimage is generated for an interventional cardiologist or a radiologist.Ideally, the DSA image mainly contains the dye-enhanced blood vesselsand they appear dark. However, the artifacts due to patient motion andsystem noise frequently reduce the diagnostic value of the images.

In the produced DSA image, a pathological aspect is associated with avascular area where there is a significant deviation from the diameterof the healthy vascular. In particular, a stenosis is associated with asignificant narrowing of the vascular and is quantified by parameterssuch as the percentage of stenosis. Stenosis limits blood flow byraising the resistance to flow through the vessel. In X-ray angiography(DSA image), the contrast agent ensures that the outlines of the bloodflow are revealed on the X-ray to indicate any narrowing of the bloodvessel. Direct visual examination of cine film coronary angiograms andmanual estimation of the degree of vascular stenosis were complicatedand subject to a large inter- and intra-observer variability.

SUMMARY

The present invention includes a system and method (i.e., utility) fordetermining the severity of a stenosis in a blood vessel. In one aspect,a method for improving DSA image quality includes: (1) registration ofthe mask and bolus images prior to subtracting procedure to reduce theartifacts from misalignment; (2) enhancement of the registered DSAimage; and (3) detecting the boundary of blood vessel and quantitativelymeasuring percentage stenosis, which may be done automatically. Thislast step is sometimes referred to as Quantitative coronary angiography(QCA). Aspects of the present invention allow such QCA to be performedby a computer with minimal user input.

The utility allows for the semi-automated quantitative measurement (QCA)of a lumen (e.g., blood vessel or artery). The utility involves aregistration step of the mask and the bolus images prior to subtractingprocedure and QCA measurement. Generally, such registration is a motioncompensation step can reduce the artifacts from misalignment and improvethe accuracy of QCA measurement. As will be appreciate, in coronaryapplications movement of the heart is nearly constant and motioncompensation significantly improves overall image quality for QCApurposes.

In one arrangement, the motion compensation is performed using analgorithm that involves an inverse-consistency constrain which impliesthat the correspondence provided by the registration in one directionmatches closely with the correspondence in the opposite direction. Thismay entail a B-spline parameterization. Other image registration methodscan also be applied in the proposed system to match the mask and bolusimages.

In one arrangement, the registration is performed hierarchically using amulti-resolution strategy in both, spatial domain and in domain of basisfunctions. The registration can be performed at ¼, ½ and full resolutionusing knot spacing of 8, 16 and 32. In addition to being faster, themulti-resolution strategy helps in improving the registration bymatching global structures at lowest resolution and then matching localvessel structures as the resolution is refined. Due to thebi-directional approach, thin blood vessels edges are more prominentthereby avoiding the local minima both in iterative nature of thinvessel edge estimation and iterative nature of correction of bolusimage. The subtraction process is applied on corrected bolus images.

Image enhancement prior to QCA measurement also improves the overall QCAmeasurement. In one arrangement the registered DSA image is furtherenhanced by background diffusion to remove noise and nonlinearnormalization for better visualization. This image enhancement increasesthe contrast between the blood vessels and the background. This resultsin much improved contrast and very crisp subtraction images, in whichthe regions of interest are easily identifiable. Enhancement of theregistered DSA image can be performed using an anisotropic diffusiontechnique and a nonlinear normalization technique. Other similar imageenhancement techniques can also be applied in the proposed utility toenhance DSA images before QCA measurement.

The utility allows a user to input various types of lumen identificationa center line method where a user identifies an approximate center ofthe lumen and an edge method where a user identifies initial edges ofthe lumen. Both methods can obtain percentage stenosis automaticallywith minimal user interaction. With the initial points selected by theuser, the QCA measurement may include one or more of the following threesub-processes: (1) initial edge detection; (2) edge refinement; and (3)lesion measurement. Initial edge detection is performed where a userinputs centerline information. In any case, the utility may refine theedges of the lumen. In one arrangement, an active contour modelalgorithm may be used to refine the edges by deforming a contour to lockonto features of interest within in an image. Other edge detection orsegmentation techniques can also be applied in the system for QCAmeasurement.

For each stenosis or lesion, the following data are calculated: (a)minimum and maximum diameter; (b) lesion length (the beginning and endof the lesion are determined automatically); (c) Stenosis reference(“normal”) diameter; and (d) percent stenosis. Other measurements canalso be included in the proposed system.

To improve the overall speed of the utility, the entire softwarearchitecture may be implemented on a GPU based framework, which makesvisualization and computation faster by up to factor of 30. Thus theentire utility may implemented substantially in real-time. Further, theutility may be used with various imaging modalities and may be used with2-D or 3-D images.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates one embodiment of a X-ray DSA system.

FIG. 2 illustrates a process flow diagram for performing an improved QCAmeasurement in accordance with various aspects of the present invention.

FIG. 3 illustrates one non-limiting method for registering images.

FIG. 4 illustrates one non-limiting method for enhancing an image.

FIG. 5 illustrates a process flow diagram for performing a QCAmeasurement.

FIGS. 6A-6C illustrate an image having a vein or artery with a lesion.

FIG. 7 illustrates a process flow diagram of an edge detection,refinement and lesion measurement process.

FIGS. 8A-8C illustrate the refinement of edge boundaries of the vein orartery of FIGS. 6A-6C.

FIG. 9 illustrates a process flow sheet for an edge detection processfor detecting edge surfaces of an object in an image.

FIG. 10 illustrates a process flow sheet for an edge refinement process.

FIG. 11 illustrates a process flow sheet for a lesion measurementprocess.

DETAILED DESCRIPTION

Reference will now be made to the accompanying drawings, which assist inillustrating the various pertinent features of the various novel aspectsof the present disclosure. Although the present invention will now bedescribed primarily in conjunction with angiography utilizing X-rayimaging, it should be expressly understood that aspects of the presentinvention may be applicable to other medical imaging applications. Forinstance, angiography may be performed using a number of differentmedical imaging modalities, including biplane X-ray/DSA, magneticresonance (MR), computed tomography (CT), ultrasound, and variouscombinations of these techniques. In this regard, the followingdescription is presented for purposes of illustration and description.Furthermore, the description is not intended to limit the invention tothe form disclosed herein. Consequently, variations and modificationscommensurate with the following teachings, and skill and knowledge ofthe relevant art, are within the scope of the present invention. Theembodiments described herein are further intended to explain known modesof practicing the invention and to enable others skilled in the art toutilize the invention in such, or other embodiments and with variousmodifications required by the particular application(s) or use(s) of thepresent invention.

FIG. 1 shows one exemplary setup for a real-time imaging procedure foruse during a contrast media/dye injection procedure. As shown, a patientis positioned on an X-ray imaging system 100 and an X-ray movie isacquired by a movie acquisition system 102. An enhanced DSA image, aswill be more fully discussed herein, is generated by an enhancementsystem 104 for output to a display 106 and/or stenosis processing systemthat is accessible to an interventional radiologist.

The projection images (e.g., CT images) are acquired at different timeinstants and consist of a movie with a series of frames before, duringand after the dye injection. The series of frames include mask imagesthat are free of contrast-enhancing dye in their field of view 108 andbolus images that contain contrast-enhancing dye in their field of view108. That is, bolus frames are images that are acquired after injecteddye has reached the field of view 108. The movie acquisition system 102is operative to detect the frames before and after dye injectionautomatically to make feasible a real-time acquisition system. Oneapproach for identifying frames before and after dye injection is tofind intensity differences between successive frames, such that a largeintensity difference is detected between the first frame after dye hasreached the field of view (FOV) and the frame acquired before it.However, the patient may undergo some motion during the imageacquisition causing such an intensity difference between even successivemask images. To avoid this, the movie acquisition system 102 may alignsuccessive frames together, such that the motion artifacts areminimized. The first image acquired after the dye has reached the FOVwill therefore cause a high intensity difference with the previous framenot containing the dye in FOV. The subtraction image or ‘DSA image’obtained by subtracting a mask frame from a bolus frame (or vice versa)will contain a near-zero value everywhere if both images belong tobackground.

Generally, the subtraction image or DSA image is obtained by computing adifference between pixel intensities of the mask image and the bolusimage. The enhancement system 104 may then enhance the contrast of thesubtraction image. Such enhancement may include rescaling theintensities of the pixels in the subtraction image and/or the removal ofnoise from the subtraction image.

The acquisition and/or enhancement systems are computerized systems thatmay run application software and computer programs which can be used tocontrol the system components, provide user interfaces, and/or providefeatures of the imaging system. The software may be originally providedon computer-readable media, such as compact disks (CDs), magnetic tape,or other mass storage medium. Alternatively, the software may bedownloaded from electronic links such as a host or vendor website. Thesoftware is installed onto a hard drive and/or electronic memory of thesystem, and is accessed and controlled by the operating system. Softwareupdates are also electronically available on mass storage media ordownloadable from the host or vendor website. The software, as providedon the computer-readable media or downloaded from electronic links,represents a computer program product usable with a programmablecomputer processor having computer-readable program code embodiedtherein. The software contains one or more programming modules,subroutines, computer links, and compilations of executable code, whichperform the functions of the imaging system. The user interacts with thesoftware via keyboard, mouse, voice recognition, and otheruser-interface devices (e.g., user I/O devices) connected to thecomputer system.

FIG. 2 shows a method for improving DSA image quality and improving theaccuracy of QCA measurements that includes three primary processes: Thefirst process includes using an averaging system 202 to generate averagemask frames 204 and average bolus frames 206 from an x-ray image ormovie. These average images are motion compensated 208. Specifically,the averaged mask frame is registered with the averaged bolus frame togenerate a registered mask frame 210, which minimizes the motionartifacts before subtraction. The registered mask frame 210 is thensubtracted 212 from the average bolus frame to generate a DSA image 214.Various aspects of this process are set forth in co-pending U.S.application Ser. No. 11/609,743 entitled “Medical Image EnhancementSystem” the entire contents of which are incorporated herein byreference. The second process uses the registered DSA image 214 asinput. This is an image enhancement process 216 that uses an anisotropicdiffusion technique to reduce noise, followed by a nonlinearnormalization method to enhance the image of the blood vessel. Theresult of this is the production of an enhanced DSA image 218. The thirdprocess uses the enhanced DSA image 218 to do a QCA measurement 220 andgenerate a QCA report or value 222. The motion compensation process isdiscussed herein relation to FIG. 3, the image enhancement process isdiscussed in relation to FIG. 4 and the QCA measurement is discussed inrelation to FIGS. 5-11.

FIG. 3 illustrates the process of registering the averaged mask framewith averaged bolus frame, to minimize the motion artifacts. In thisembodiment, this is done using an inverse consistent image registrationusing B-spline basis. In medical imaging, image registration is used tofind a point-wise correspondence between a pair or group of similaranatomical objects. Image registration is often posed as an optimizationproblem that minimizes an objective function representing the differencebetween two images 302, 304 to be registered. In this case, the imagesare the average bolus and mask frames. The symmetric squared intensitydifference is chosen as the driving function. In addition,regularization constraints are applied so that the deformation follows amodel that matches closely with the deformation of real-world objects.The regularization is applied in the form of bending energy andinverse-consistency cost. Inverse-consistency implies that thecorrespondence provided by the registration in one direction matchesclosely with the correspondence in the opposite direction. Most imageregistration methods are uni-directional and therefore containcorrespondence ambiguities originating from choice of direction ofregistration. In the presented method, the forward and reversecorrespondences are computed simultaneously and bound together throughan inverse consistency cost term. The inverse consistency cost termassigns higher cost to transformations deviating from beinginverse-consistent. While inverse consistency minimizes thecorrespondence ambiguity, it also helps the transformation performbetter by forcing it out of local minima. A cost function for performingimage registration over the images is used:

$\begin{matrix}{C = {{\sigma \begin{pmatrix}{{\int_{\Omega}{{{{I_{1}\left( {h_{1,2}(x)} \right)} - {I_{2}(x)}}}^{2}{x}}} +} \\{\int_{\Omega}{{{{I_{2}\left( {h_{2,1}(x)} \right)} - {I_{1}(x)}}}^{2}{x}}}\end{pmatrix}} + {\rho \begin{pmatrix}{{\int_{\Omega}{{{L\left( {u_{1,2}(x)} \right)}}^{2}{x}}} +} \\{\int_{\Omega}{{{L\left( {u_{2,1}(x)} \right.}^{2}{x}}}}\end{pmatrix}} + {\chi \begin{pmatrix}{{\int_{\Omega}{{{{h_{1,2}(x)} - {h_{2,1}^{- 1}(x)}}}^{2}{x}}} +} \\{\int_{\Omega}{{{{h_{2,1}(x)} - {h_{1,2}^{- 1}(x)}}}^{2}{x}}}\end{pmatrix}}}} & (1)\end{matrix}$

where, I₁(x) and I₂(x) represent the intensity of image at location x, xrepresents the domain of the image. h_(i,j)(x)=x+u_(i,j)(x) representsthe transformation that maps image I_(i) to image I_(j) in the Eulerianframe of reference and u(x) represents the displacement field. L is adifferential operator and the second term in Eq. (1) represents anenergy function. σ, ρ and χ are weights to adjust relative importance ofthe cost function.

In Equation (1), the first term represents the symmetric squaredintensity cost function and represents the integration of squaredintensity difference between deformed reference image and the targetimage in both directions. The second term represents the energyregularization cost term and penalizes high derivatives of u(x). L isrepresented as a Laplacian operator mathematically given as: L=∇². Thethird and last term represents the inverse consistency cost function,which penalizes differences between transformation in one direction andinverse of transformation in opposite direction. The total cost 308 iscomputed 306 using Eq. 1 as a first step in registration.

The optimization problem posed in Eq. (2) is solved by using a B-splineparameterization as is known in the art. B-splines are used due to easeof computation, good approximation properties and their local support.It is also easier to incorporate landmarks in the cost term if we usespatial basis function. The above optimization problem is solved bysolving for B-spline coefficients c_(i)'s, such that

$\begin{matrix}{{h(x)} = {x + {\sum\limits_{i}{c_{i}{\beta_{i\;}(x)}}}}} & (2)\end{matrix}$

where, β_(i)(x) represents the value of B-spline at location x,originating at index i. In this registration method, cubic B-splines areused. A gradient descent scheme is implemented based on the aboveparameterization. The total gradient cost is calculated and updated 310with respect to the transformation parameters in every iteration. Thetransformation parameters are updated using the gradient descent updaterule. Images are deformed into shape of one another using the updatedcorrespondence and the cost function and gradient costs are calculated314 until convergence 316 when the frames are registered 318. Theregistration is performed hierarchically using a multi-resolutionstrategy in both, spatial domain and in domain of basis functions. Theregistration is performed at ¼, ½ and full resolution using knot spacingof 8, 16 and 32. In addition to being faster, the multi-resolutionstrategy helps in improving the registration by matching globalstructures at lowest resolution and their matching local structures asthe resolution is refined.

Referring again to FIG. 2, these registered frames are subtracted toprovide the DSA image for image enhancement. That is, after imageregistration, subtraction is performed between the registered mask imageand bolus image, resulting in the motion compensated (registered) DSAimage 214. Using the registered DSA image as input, the system furtherimproves image quality by the following two steps: 1) noise reductionthrough an anisotropic diffusion technique that can effectively removethe background noise without destroying the vessels and translucency ofthe contents inside the vessels; and 2) enhance the blood vessel imagethere by suppressing the background through nonlinear normalization.This procedure is shown in FIG. 4.

Initially, the DSA image 214 is provided for look-up-table (LUT) baseddiffusion 402 or anisotropic diffusion. Such diffusion 402 is describedin relation to equations 3-6. Such an anisotropic diffusion and is basedon partial differential equation (PDE) for noise smoothing. Given animage l(x,y,t) at time scale t, the diffusion process is expressed as:

$\begin{matrix}{{\frac{\partial}{\partial t}{I\left( {x,y,t} \right)}} = {{div}\left( {{c\left( {x,y,t} \right)}{\nabla I}} \right)}} & (3)\end{matrix}$

where ∇ is the gradient operator, div is the divergence operator, andc(x,y,t) is the diffusion coefficient at location (x,y) at time t, Withapplying the divergence operator, Eq. (3) can be rewritten as

$\begin{matrix}{{\frac{\partial}{\partial t}{I\left( {x,y,t} \right)}} = {{{c\left( {x,y,t} \right)}\Delta \; I} + {{\nabla{c\left( {x,y,t} \right)}}{\nabla I}}}} & (4)\end{matrix}$

where Δ is the Laplacian operator. The diffusion coefficient c(x,y,t) isthe key in the smoothing process and it should encouragehomogenous-region smoothing and inhibit the smoothing across theboundaries. It is chosen as a function of the magnitude of the gradientof the brightness function, i.e.

c(x,y,t)=g(∥∇I(x,y,t)∥)  (5)

The suggested functions for g(·) are the following two

$\begin{matrix}{{g\left( {\nabla I} \right)} = {{^{- {(\frac{{\nabla I}}{K})}^{2}}\mspace{20mu} {or}\mspace{14mu} {g\left( {\nabla I} \right)}} = \frac{1}{1 + \left( \frac{{\nabla I}}{K} \right)^{2}}}} & (6)\end{matrix}$

where K is the diffusion constant which controls the edge magnitudethreshold. Generally speaking, a larger K produces a smoother result ina homogenous region than a smaller one. Here we apply diffusiontechnique on the input DSA images to smooth background thereby reducingthe structured noise. Also, the algorithm may be implemented usingcompute unified device architecture (CUDA), which is developed by nVidiato run general purpose computations on a Graphics Processing Unit (GPU).The result of the diffusion process 402 is a diffused image 404.

To increase the contrast between the blood vessels of and background inthe DSA image, a nonlinear normalization method 406 is used. Theprinciple is to force the background to suppress the non-vascularstructures and noise, while gradually enhancing the foreground vascularstructures. The nonlinear normalization 406 is currently implemented asa LUT filter that is based on the pre-defined parameters and is setforth in relation to equation 7.

Letting I_(in)(x,y) be the input DSA image (after diffusion), thenonlinear normalization, defined as I_(out)(X,Y) is mathematically givenas:

$\begin{matrix}{{I_{out}\left( {x,y} \right)} = \left\{ \begin{matrix}{{I_{t}\left( \frac{I_{i\; n}\left( {x,y} \right)}{I_{t}} \right)}^{y_{1}},{{I_{i\; n}\left( {x,y} \right)} \in \left\lbrack {0,I_{t}} \right\rbrack}} \\{{\left( {I_{t} + 1} \right) + {\left( {I_{t} + 1} \right)\left( \frac{{I_{i\; n}\left( {x,y} \right)} - \left( {I_{t} + 1} \right)}{I_{t} + 1} \right)^{y_{2}}}},{{I_{i\; n}\left( {x,y} \right)} \in \left\lbrack {I_{t} + {1\text{,}255}} \right\rbrack}}\end{matrix} \right.} & (7)\end{matrix}$

Here I_(i) is a pre-defined threshold, y₁>1.0 and y₂<1.0 for class-basedcontrast enhancement. So if the intensity range of dye lies in lowerhalf of the image and the background lies in the higher half of theimage, the blood vessels are enhanced in limits. The result of thisnon-linear normalization is the enhanced DSA image 218.

After these processes, the DSA image quality has been improved greatly.Such improvement in the image facilitates QCA measurement as set forthherein. However, it will be appreciated that the QCA measurement setforth herein is considered novel in and of itself.

Set forth herein are two semi-automatic algorithms to compute percentagestenosis with minimum user intervention. As shown in FIG. 5, the overallprocess flow for QCA is as follows. Using a DSA image, which may be anenhanced image as set forth above, a user or operator may select 502 oneframe from the DSA to use in stenosis estimation. This beings the QCAprocess 504. Once the frame is selected a user selects a region ofinterest 506. This may entail zooming in 508 on the region of interestto provide an enlarged view of the region of interest 510. This isillustrated in FIG. 6A which shows a frame having a blood vessel 600that includes a lesion region 602 or stenosis.

The operator can then choose “measure from centerline” (center linemethod 512), or “measure from sides” (vessel edge method 514). If“measure from centerline” is chosen the operator needs to select 516several points 604 along the centerline of vessel 600 as illustrated inFIG. 6B. If “measure from sides”, is chosen the operator needs to select516 several points 604 along the both sides of vessel 600 as illustratedin FIG. 6C. The QCA process 518 is then performed by the system whichsubsequently generates a final report 222 that indicates a percentage ofstenosis. The final stenosis results is computed and displayedautomatically.

The QCA measurement process is set forth in FIG. 7. With the initialpoints selected by the user, the QCA measurement includes the followingthree sub-processes: (1) initial edge detection 720; (2) edge refinement740; and (3) lesion measurement 760.

When the user chooses initial points 604 along the centerline 610,initial edge detection in QCA is performed by calculating the gradientsof in a direction perpendicular to the centerline across the lumen andfinding their peak values. See FIG. 8A. As will be appreciated, thecontrast between the dyed vessel and edge is usually significant andeasily detectable. That is, the vessel is often dark or nearly black inan image whereas the vessel wall and surrounding area is grey or white.This process is performed along the centerline 610 until edge points 612are determined along both edges of the vessel 600 through the stenosisregion 602. See FIG. 8B

FIG. 9 shows the process of initial edge detection. Initially, after thecenterline points 604 are selected and connected to form the centerline610, perpendicular line construction 902 is performed to form linesperpendicular 904 to the centerline. On each perpendicular line, theweighted sum of the first and second gradients is calculated 906 for thesubject target line. The peak value is used to find the resultinginitial boundary or initial edge points 612. These initial edge pointsare detected by using Eq. (8).

$\begin{matrix}{{f_{edge}\left( {x,y} \right)} = {\underset{f{({x,y})}}{\arg \mspace{11mu} \max}\left( {{w_{1}{\nabla{f\left( {x,y} \right)}}} + {w_{2}{\nabla^{2}{f\left( {x,y} \right)}}}} \right)}} & (8)\end{matrix}$

If the user chooses initial points along the two boundary lines, thosepoints are used as the initial edge points and above initial edgedetection step is omitted.

After initial edge point detection, a group of edge points areidentified along two boundaries of the vessel. See FIG. 8B. Theseinitial points are refined 740 (see FIG. 7) to get a smooth boundarycurve. This process is set forth in FIG. 10. This process uses an activecontour model algorithm, which will be known to one skilled in the art,that deforms a contour to lock onto features of interest within in animage. Usually the features are lines, edges, and/or object boundaries.Given an approximation of the boundary of an object in an image (e.g.,initial point 612), an active contour model can be used to find the“actual” boundary. An active contour is an ordered collection of npoints in the image plane:

V={v₁, v₂ . . . v_(n)}

v _(i)=(x _(i) ,y _(i)), i=1, . . . n  (9)

The points in the contour iteratively approach the boundary of an objectthrough the solution of an energy minimization problem. For each pointin the neighborhood of v_(i), an energy term 902 is computed:

E _(i) =αE _(int)(v _(i))+βE _(ext)(v _(i))  (10)

where E_(int)(v_(i)) is an energy function dependent on the shape of thecontour and E_(ext)(v_(i)) is an energy function dependent on the imageproperties, such as the gradient, near point v_(i). α and β areconstants providing the relative weighting of the energy terms.

The internal energy function is intended to enforce a shape on thedeformable contour and to maintain a constant distance between thepoints in the contour. Additional terms can be added to influence themotion of the contour. The internal energy function used herein isdefined as follows:

αE _(int)(v _(i))=bE _(con)(v _(i))+cE _(cur)(v _(i))  (11)

where E_(con)(v_(i)) is the continuity energy that enforces the shape ofthe contour and E_(cur)(v_(i)) is a curvature energy that causes thecontour to grow or shrink. c and b provide the relative weighting of theenergy terms.

The external energy function attracts the deformable contour tointeresting features, such as object boundaries, in an image. Here imagegradient is used. The image gradient should be large at the objectboundary (β, <0). Therefore, the following external energy function isinvestigated:

βE _(int)(v _(i))=βE _(grad)(v _(i))  (12)

In summary, energy function at each v_(i) is minimized:

E _(i) =bE _(con)(v _(i))+cE _(cur)(v _(i))+βE _(grad)(v _(i))  (13)

In (8), b>0,c>0 and β<0.

This process iteratively updates the initial edge points 904 until theenergy function is minimized or a maximum number of iterations areachieved The result are refined edge points 910 as illustrated in FIG.8C. That is, the points are smoothed and can be fitted with a curve.

After these first two processes 720, 740, the edges of the lumen areidentifies a dimensions need to be calculated for lesion measurement740. See FIG. 7. If dealing with a lesion, it is usually desirable toidentify at least the minimum luminal diameter (MLD), the percentstenosis and perhaps the length of the lesion.

One method to determine lumen dimensions is set forth in the processflow sheet of FIG. 11. Using the detected boundaries 1102 as set forthby the refined edge points, a centerline is constructed 1104 between thetwo edges that represents the initial centerline 1106 of the lumen.Diameters are then measured by constructing lines that are perpendicularthe centerline. Then, it is determined where each diameter lineintersects 1108 each of the two edges. This defines a set ofcorresponding points on the boundary 1110. The distance is calculated1112 between the two points of intersection for each line. The distancehaving the shortest length is the minimum luminal diameter. The MLD maybe divided by the normal diameter of the vessel to provide a stenosispercentage 1114. Generally, for each lesion, the following data arecalculated: (a) Minimum and maximum diameter; (b) Lesion length (thebeginning and end of the lesion are determined automatically); (c)Stenosis reference (“normal”) diameter; and (d) Percent stenosis. In anycase, this information may be output to the operator.

The proposed system provides a number of advantages. For instance, maybe used in the quantitative measurement of a blood vessel and the systemincludes three sub-systems: motion compensation, image enhancement, andstenosis measurement. The whole system results in a more intelligent andmore accurate system for improving QCA measurement. Another advantage isthat the proposed system can involve a registration step to align themask and bolus images prior to subtracting procedure and QCAmeasurement. This motion compensation step can reduce the artifacts frommisalignment (due to patient movement) and improve the accuracy of QCAmeasurement. Another advantage is that the proposed system can involvean image enhancement step prior to QCA measurement; the resulted DSAimage is further enhanced by background diffusion and nonlinearnormalization for better visualization. This image enhancement increasesthe contrast between the blood vessels and the background. This resultsin much improved contrast and very crisp subtraction images, in whichthe regions of interest are easily identifiable.

It will be appreciated that the entire software architecture in theproposed system may be implemented on a GPU based framework, which makesvisualization and computation faster by up to factor of 30. The entirescheme may therefore be implemented as a real-time scheme.

The foregoing description of the present invention has been presentedfor purposes of illustration and description. Furthermore, thedescription is not intended to limit the invention to the form disclosedherein. Consequently, variations and modifications commensurate with theabove teachings, and skill and knowledge of the relevant art, are withinthe scope of the present invention. The embodiments describedhereinabove are further intended to explain best modes known ofpracticing the invention and to enable others skilled in the art toutilize the invention in such, or other embodiments and with variousmodifications required by the particular application(s) or use(s) of thepresent invention. It is intended that the appended claims be construedto include alternative embodiments to the extent permitted by the priorart.

1. A method for use in the quantitative measurement of a blood vesselblockage, comprising: operating an imaging system to obtain a series ofmask images and a series of bolus images corresponding to a region ofinterest; operating a processor to register the mask images to the bolusimages to generate a registered mask image and subtract the registeredmask image form the bolus image to produce a motion compensated DSAimage; operating the processor to process the DSA image to enhancecontrast in at least the region of interest of the DSA image; receivinga user input identifying a lumen within the region of interest; andbased on the user input, operating the processor to autonomouslyidentify a boundary of the lumen and calculate a stenosis percentage forthe lumen.
 2. The method of claim 1, wherein registering the series ofmask images and the series of bolus images further comprises averagingsaid series to produce an average mask image and an average bolus image,wherein the average mask image is registered to the average bolus image.3. The method of claim 1, wherein registering comprises applying aninverse-consistency constraint to a deformation of the mask frame to thebolus frame.
 4. The method of claim 3, wherein the inverse-consistencyconstraint comprises a B-spline parameterization.
 5. The method of claim1, wherein enhancing the region of inertest comprises: operating theprocessor to diffuse the DSA image to remove background noise from theDSA image.
 6. The method of claim 5, further comprising: operating theprocessor to perform a non-linear normalization on the DSA image,wherein a contrast between the lumen and the background of the DSA imageis increased.
 7. The method of claim 1, wherein receiving the user inputidentifying the lumen within the region of interest comprises one of:receiving at least two user selected points associated with a centerlineof said lumen; receiving a plurality of user selected points associatedwith edges of the lumen.
 8. The method of claim 7, wherein operating theprocessor to autonomously identify a boundary of the lumen and calculatea stenosis percentage for the lumen further comprises at least one of:initial edge detection; edge refinement; and lesion measurement.
 9. Themethod of claim 8, wherein initial edge detection based on the at leasttwo points associated with the centerline of the lumen, comprisesdefining an initial centerline in the lumen; projecting a plurality ofperpendicular lines relative to the centerline; and for eachperpendicular line, identifying a boundary point between the lumen andthe background based on gradients calculated along said perpendicularline, wherein a plurality of the boundary points define an initial edgeof the lumen.
 10. The method of claim 8, wherein edge refinementcomprises fitting a contour to the initial edge, wherein the contour issmoothed to define a refined edge.
 11. The method of claim 8, whereinlesion measurement comprises: using the processor to calculate (a)minimum and maximum diameters of said lumen; (b) lesion length; (c) astenosis reference (“normal”) diameter; and (d) percent stenosis. 12.The method of claim 1, wherein the steps of operating a processorcomprise operating a GPU based processor.
 13. The method of claim 1,wherein said series of mask and bolus frames are acquired using at leastone of the following imagining modalities: X-ray; CT; MRI; andultrasound.
 14. The method of claim 1, wherein operating a processor toregister the mask image to the bolus images to generate a registeredmask image comprises: registering a 3-D mask image to a 3-D bolus imageto generate a 3-D registered mask image.
 15. The method of claim 1,wherein operating a processor to register the mask image to the bolusimages to generate a registered mask image comprises: performing aninitial registration at a reduced resolution to align global structures;and performing a secondary registration to align local structures in theregion of interest.
 16. A system for use in the quantitative measurementof a blood vessel blockage, comprising: an imaging system operative toobtain a series of mask images and a series of bolus imagescorresponding to a region of interest; an image processing systemoperative to: register the mask image to the bolus images and generate aregistered mask image; subtract the registered mask image form the bolusimage to produce a motion compensated DSA image; process the DSA imageto enhance contrast in at least the region of interest of the DSA image;receive a user input identifying a lumen within the region of interest;and based on the user input, operating the processor to autonomouslyidentify a boundary of the lumen and calculate a stenosis percentage forthe lumen.
 17. The system of claim 16, wherein the imaging systemcomprises an X-ray imaging system.
 18. A method for use in thequantitative measurement of a blood vessel blockage, comprising:operating an imaging system to obtain a series of mask images and aseries of bolus images corresponding to a region of interest; operatinga processor to register the mask images to the bolus images to generatea registered mask image and subtract the registered mask image form thebolus image to produce a motion compensated DSA image; receiving a userinput identifying a lumen within the region of interest; and based onthe user input, operating the processor to autonomously calculate: (a)minimum and maximum diameters of said lumen; (b) lesion length; (c) astenosis reference (“normal”) diameter; and (d) percent stenosis of thelumen.
 19. The method of claim 18, further comprising identifying edgesfor the lumen by: defining an initial centerline in the lumen;projecting a plurality of perpendicular lines relative to thecenterline; and for each perpendicular line, identifying a boundarypoint between the lumen and the background based on gradients along saidperpendicular line, wherein a plurality of the boundary points define anedge of the lumen.